Low-cost fermions in classical field simulations 
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We discuss the possible extension of the bosonic classical field theory simulations to include 
fermions. This problem has been addressed in terms of the inhomogeneous mean field approximation 
by Aarts and Smit. By performing a stochastic integration of an equivalent set of equations we can 
extend the original 1+1 dimensional calculations so that they become feasible in higher dimensions. 
We test the scheme in 2 + 1 dimensions and discuss some classical applications with fermions for 
the first time, such as the decay of oscillons. 
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I. INTRODUCTION 

Since the advent of modern computational facilities 
classical field theory is one of the most popular ap- 
proaches to nonequilibrium field theory. The classical 
approximation to a quantum filed theory is well justi- 
fied in several cosmological applications ranging form re- 
heating of the postinfiationary Universe [l|, 01 followed 
by an evolution through various phase transitions [3| to 
the nonlinear evalution of the hypothetical cosmic strings 
01 S IB'] . Classical methods have also received an increas- 
ing amount of attention from the heavy ion community. 
The initial evolution of the highly excited gluon plasma 
in little-bang experiments turns out to be well modelled 
by classical Yang-Mills equations ^Ti] . 

The preheating of the inflationary Universe was one 
of the pioneering applications of the nonlinear classical 
wave equations [8j] . In the mostly studied chaotic and hy- 
brid inflation scenarios the nonlinear dynamics is driven 
by an instability, which is parametric or tachyonic, re- 
spectively. Instabilities lead to nonpcrturbatively large 
occupation numbers, which is a prerequisite for the clas- 
sical approximation, but it also requires a nonperturba- 
tive treatment, which is the actual strength of the classi- 
cal equations. The classical simulations of preheating can 
make estimates on non-gaussian density perturbations ^] 
and gravitational wave production [l^, [ill [I4I . 

Nonperturbative methods are especially useful when 
dealing with nonequilibrium phase transitions in the 
Early Universe [l3]. A typical example where the fields 
are required to be out of equilibrium is baryogenesis. A 
second strength of the classical approximation is that 
equilibrium is not a prerequisite. Solving the real time 
Yang-Mills equations with far-from-equilibrium initial 
conditions one could gain access to the evolution of the 
Chern-Simons number , llBl ■ For this to accomplish 

the third strength of the classial equations has been ex- 
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ploited: its preservation of gauge invariance under time- 
independent transformations. 

The fourth strength of the classical approach is its sim- 
plicity and cheap implementation even at large scales. 
This feature makes it an excellent tool for studying topo- 
logical defects, especially the hypothetical network of cos- 
mic strings. To address formation and evolution of defect 
networks very different length scales have to be properly 
incorporated into one numerical computation. This sit- 
uation is getting worse in an expanding universe, but in 
the classical setting these calculations are still affordable 
In principle, one could take the zero-width limit 
and solve the Nambu-Goto equations For funda- 

mental strings this is a natural procedure, but for strings 
which are topological defects, microscopic physics plays 
a significant role in the decay mechanism of strings [l9|. 
Explicit calculations have been made in the context of 
gauge strings in the Abelian Higgs model [13, [3 l20l |. 
global [2l| and semilocal strings [2^ . [23l | as well as do- 
main walls [13, l25J . The relevance of these calculations 
have been recently highlighted by the discovery of pos- 
sible traces of cosmic strings in the cosmic microwave 
background [23l. [26l|. 

Besides of the cosmological interest classical field the- 
ory simulations are also used on the subatomic scale: 
the early evolution of the gluon plasma formed in heavy 
ion collisions can be described by by classical Yang-Mills 
equations [3|- This facilitates a nonperturbative descrip- 
tion of the glasma, i.e. the intermediate state after the 
melting of the color glass condensa te p rior to thermal- 
isation to quark gluon plasma [13, [13 . The produced 
nonabelian plasma is highly anisotropic, and as such, it 
is subject to instabilities f29'|. Classical methods have 
proved very useful for giving a quantitative account on 
the isotropisation driven by these Weibel instabilities 
[30I [sTj . Alternatively, one can replace the hard sector of 
the field theory by classical particles represented by a set 
of Vlasov equations on the background of soft classical 
fields [3^. [33^. 

Finally we point out that the ergodicity of the classical 
field trajectories makes the classical simulations an es- 
sential and robust method for studying thermal classical 
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lattice systems in real time. In statistical field theory one 
averages over an ensemble of initial field configurations 
and observes e.g. the real time dynamics of symmetry 
breaking [sJl, with possible formation of quasistable lo- 
calised excitations, dubbed oscillons [35J. The presence 
of long-lived oscillons induces resonant nucleation, and 
they become a driving force of first order phase transi- 
tions [36]. 

The classical approximation has severe limitations, 
however. The continuum equilibrium theory is plagued 
by Rayleigh- Jeans divergences, and a renormalisation 
with local counterterms is not possible in general [svj . 
Moreover, the counterterms are temperature dependent, 
which makes a consistent out-of-equilibrium renormalisa- 
tion impossible. This means that classical theories need 
an intrinsic cut-off scale, which, in practice, sets the spac- 
ing of the lattice discretisation. From whatever initial en- 
semble of classical fields the straightforward integration 
of the Euler-Lagrange equations of the theory brings the 
systems towards an equilibrium defined by the classical 
Hamiltonian. This equilibrium differs from a true quan- 
tum thermal state, but the difference is negligible for 
soft modes and only affects hard excitations. In terms 
of particle numbers a system is considered in the clas- 
sical domain if the occupancy is sufficiently high. The 
infrared physics, which is mostly sensitive to nonpertur- 
bative phenomena, is usually not vulnerable to quantum 
effects, but on the ulraviolet end of the spectrum one 
has to balance between discretisation errors and miscal- 
culated hard degrees of freedom. Even if we start from 
an infrared dominated initial condition, hard modes are 
becoming increasingly dominant on the course of ther- 
malisation and the classical system automatically leaves 
its domain of validity. 

There is an other first principles approach to nonequi- 
librium field theory, which shares none of the aforemen- 
tioned shortcomings. It has been numerically demon- 
strated that even a low order truncation of the two- 
particle irreducible (2PI) effective action yields equations 
of motion, capable of describing irreversible quantum dy- 
namics, including thermalisation (ssj . This powerful re- 
summation technique can be directly applied to relevant 
problems in cosmology [H, [13] or in hot abelian gauge 
theories [iH , as well as in the many-body theory of 
ultracold condensates [i^. Yet, for non- abelian gauge 
fields the more complete 3PI resummation becomes nec- 
essary m, |45|] , and for a setting with topological defects 
an inhomogeneous treatment is inevitable |46|. Being 
both extensions expensive, we will have to fall back in 
these cases to the classical approximation and use 2PI 
to benchmark it where their domains of validity overlap. 
These precision tests in the 0{N) scalar model had the 
reassuring result that particle numbers as small as ~ 10 
already put the system into the classical domain [i^, . 

There is, however, another important deficiency of the 
classical approximation. All the applications listed in 
the previous paragraphs were entirely limited to bosonic 
fields. Classical simulations of both baryogenesis and 



heavy ion collisions could benefit from a direct modelling 
of quarks, if that were feasable. 

As dimensional reduction suggests, fermionic fields are 
purely quantum degrees of freedom, just like the non- 
static components of a bosonic field theory. The classical 
field theory does have bosonic fluctuations, and in the 
absence of quantum degrees of freedom, bosonic particle 
production is automatically modelled as the excitation of 
the fluctuating background. The analagous production 
of fermions, however, is not mapped onto any existing 
degree of freedom. 

The inclusion of fermions is rather trivial in the 2PI 
framework, where bosonic quantum fluctiatons interact 
with fermionic quantum fluctuations, and an explicit cal- 
culation has been presented to show the real-time simul- 
taneous onset of Fermi-Dirac and Bose-Einstein distri- 
butions [4^ . However, when dealing with non-abelian 
gauge fields, or strong inhomogeneities, we will need to 
resort to some extension of the classical theory. This 
extension is the actual topic of this paper. 

In this paper we build on the ideas of Aarts and 
Smit [13, [5l| and by "integrating the fermion determi- 
nant" we solve an effective theory for the classical scalar 
background. We go beyond the recent applications in 
Refs. m, HI] by including the back reaction in our cal- 
culation. Our efficient solution technique enables us to 
go beyond 1-1-1 dimensions in the simulations. 

In Section |ll] we review the standard description of 
the fermionic fluctuations. Then in Section IIIII we in- 
troduce a stochastic approach, which provides us a more 
efficient algorithm than the so far known mode function 
expansion. In Section HVl wc investigate the capabilities 
of this semiclassical approximation for describing irre- 
versible phenomena, such as damping and thermalisa- 
tion. We continue with a bit more exotic application 
involving oscillons in Section |V] and discuss the possible 
future applications of this semiclassical scheme Section 
IVII The spinor representations that we actually used in 
our numerics we give in Appendix [Xj In a naively discer- 
tised lattice field theory the number of fermion flavours 
is doubled in each space-time direction. We discuss the 
possible elimination of the extra flavours in Appendix IB] 



II. INTEGRATING THE FERMION DEGREE 
OF FREEDOM 

A. A scalar model with fermions 

Let us pick a simple scalar model coupled to a fermion 
flavour through Yukawa interaction: 

+ E b^^hrd^^k - ^k{MPL + M*Pn)^k] (1) 

k 
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Here the M complex fermion mass is a function of the 
background: 

M{x) = m- g^{x) . (2) 

The projectors are defined as Pl = ^(1 — 7^) and Pu = 
^{1 + 7^). The index k runs over Nf identical fermion 
flavours. We will not use any of the special features of 
the bosonic sector, and our discussion below will also 
apply to classical lattice gauge theories with a covariant 
coupling to fermions. 

Before going into details we summarise our strategy by 
defining a bosonic effective action r[$] as 

^^^m ^ I J]i5vi/+DvI/fee^/^(*'*+'*) . (3) 

k 

Our goal is to solve the semiclassical equation of motion 
(5r[$]/(5<f> — without further approximation. This path 
integral has to be understood on a real time contour with 
a forward and backward time branch. To contour ends 
at (zero) initial time where it connects to the initial den- 
sity operator. We will use a the perturbative vacuum or 
many-particle state as an initial condition. 

The Dirac equation written for the spinor operators is 
linear 

{il'^df, -711 + gRe^[x) - ?s(Im$(x)7^)*(x) = 0,(4) 
«a^*(a;)7^ + ^{x)[m - gRe$(x) + i.glm $(a;)7^) = 0,(5) 

which manifests on the level of diagrams in the simple 
rule that fermion propagator lines never cross. The inter- 
action is mediated by the bosonic field, which is modelled 
by a fluctuating background. 

Instead of using anticommuting operators we rewrite 
the Dirac equation so that it applies to the symmetrised 
two-point function: 

{ij'^dx,!! - rn 

+gRe ^{x) - igliR ^{x)j^)D{x, y) 0, (6) 
idy,^D{x,y)-f^' + 

D{x,y)im^ gRe $(y) + iglui $(y)7^) = , (7) 

where D(x,y) is defined as 

D{x,y),, = i(D>(x,2/)-Z?<(x,y)) 

= i([vl/,(x),^,(y)]) , (8) 

D<(x,y\, = (*,(y)vl/,(x)) , (9) 
D>{x,y\, - {^.,{x)^,{y)) . (10) 

where i,j represent the Dirac as well as flavour indices. 
The propagator D is identical to the F-type two-point 
function in the literature of nonequilibrium Green's func- 
tions as well as in Ref. 4!!]. One can work out an equa- 
tion for the spectral function as well, which will take an 
identical form. 



The bosonic background obeys a simple wave equation, 

dl<^ix) +V'{^{x)) + NfJ{x) =0 (11) 

where the fermionic back reaction is carried by the cur- 
rent J, which is a combination of the scalar and pseu- 
doscalar currents: 

J{x) ^ J^{x) + J^^{x) = 2gTTD{x,x)PR , (12) 
J^{x) = ~g{^{x)-i'{x)) = gTTD{x,x) , (13) 
jP^(x) = -g{^'{x)-/^^{x)) ^ gTTD{x,x)-f^ . (U) 

The scalar current is always real the pseudoscalar current 
is always imaginary. 

In a theory with a Dirac mass m the vacuum propaga- 
tor takes the following form 

D{x',xy,y)U^,-l^e'^''^'''-y'^^^^, (15) 

with uj'i = + The latin indices refer to space 

only. In this equation we introduced the notation J- for 
the three dimensional momentum integral J dPp/ {2Tr)^ . 

In many cases when one inquires about the fermion 
production the vacuum initial condition is used, prefer- 
ably. Since Eqs. are flrst order in time, all fur- 
ther evolution is determined, once the background is 
known. Of course, any other initial particle content is 
also feasable, one can e.g. set an uneven number of par- 
ticles and antiparticles, which is the microcanonical ana- 
log of a baryochemical potential. We will give formulas 
where these particle numbers enter later below. 

B. Mode function expansion 

One can solve Eqs. ([B][7]) and pT|l numerically without 
any further information. The standard strategy is to in- 
troduce mode functions, i.e. to treat time evolution as a 
Bogolyubov transformaton of the initial-time ladder op- 
erator. This method has been formerly used for bosonic 
fluctuations on a homogeneous background [sS , [ssl . [5^ , 
and later extended to fermionic systems [stI . l58j and also 
to inhomogeneous backgrounds [HI, [gO] ■ The equations 
for fermionic fluctuations on an inhomogeneous back- 
grounds have been worked out in detail by Aarts and 
Smit [13 . 

We introduce the mode functions (j)^'''{x,p) and 
(f>^''^{x,p) as classical solutions weighting the anticom- 
muting ladder operators with 

{bM,bt'iQ)} = {^^f5{j>~ q)5,,s' , (16) 
{dM,d+{q)] = {27rf S{p~ q)S,,,, (17) 

in the fermion held operator: 

(18) 
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We introduced the spinor index s that runs from 1 to 2. 
If the fcrmions' initial condition is homogeneous, one has 



(19) 
(20) 



The ladder operators correspond to these initial time ex- 
citations that are transformed as fermions travel through 
the background. The statistical features of these opera- 
tors actually reflect the initial particle distribution: 



These equations relate the propagators used in Eqs. ([6][7|) 
to the mode functions. So that '^{x) in Eq. solves 
the Dirac equation ^ the mode functions (x,7?) as 
well as (f)^''^ {x,p) have to solve the same Dirac equation 
for all p and s. 

{ij^df, -m + gRe$(x) i5lm$(2;)7^ )</)"/"'" (x,p) =0, 

(33) 

This is now also manifest from Eq. ([5^ . We can actu- 
ally confirm the initial condition in Eq. (|20|) by inserting 
D{x,y) of Eq. ^ into Eq. 



¥{p),b^ (q) ) - {27rf S{p- q)Ss,s'a ^ 2<(p)) , (21) 



d'{p),d^ {q) ) = {2nY5{p- q)5s^s'{l - 2711 (p)) . (22) 



The u'^{p) and v'^{p) spinors in Eq. (|20)) are defined as 
the eigenvectors of the vacuum correlation matrix written 
momentum space: 



1 



(23) 



This matrix has the eigenvalues (+1, +1, — 1, — 1) corre- 
sponding to the eigenvectors u^{p), u^ip), v^{p) and v'^{p) 
respectively. On a non-trivial background these eigenval- 
ues disambiguate between particle and antiparticle solu- 
tions. Using the identities 



C. Renormalisation 

The effective potential in Eq. ([3]) has a non-polynomial 
contribution from the logarithm of the fermion determi- 
nant. Expanding in $ to n-th order one finds the fermion 
one-loop diagrams with n external bosonic lines. These 
diagrams with n < 4 are potentially divergent in 3-1-1 di- 
mensions. Already at ti = 1, the source is quadrati- 
cally divergent. 

We renormalise the scalar potential additively by intro- 
ducing a renormalised potential V and a counterfunction 
SV'{^) in Eq. (jlip . We also introduce a wave function 
renormalisation so that the renormalised scalar evolution 
equation reads 



= ys^^ ^ (24) Zdl^Rix) + Vl,i<!>R{x)) + SV'i^Rix)) + NfJ{x) = . 



S 

J2 {u%p)u^+{p) - v%p)v^ + ip)) = Mp (26) 



one can show that at initial time the two-point function 
in Eq. (fTS]) is correctly reproduced by the field operator 
in Eq. (UHl). 

At any later xq the mode functions are given by the 
following commutators: 

{[^{x),b^+ip}]) = r'^(x,p), (27) 
{[^ix),d^^p)]) = -0^'^(x,p). (28) 

On the other hand, one can express the ladder operators 
in terms of the initial time field operator by 



(34) 

To calculate SZ = Z — 1 we linearise J in $ and obtain 

Zdl^Rix) + VI,{^r{x)) + 6r{<S>R{x)) = 

Nf J dzo j d^zY.{x~ z)^r{z) . (35) 

Here S(a;) stands for the vacuum one-loop self energy. 
The counterterms 6Z and (see Eq. PO)) below) will 
be set so that they cancel the potentially divergent first 
two coefficients in the k'^ expansion of Yi{kQ,k) so that 
the renormalised self energy 



Y.R = SZk^ + 5m^ + E(fco, k) 



(36) 



¥^{p) = / vI/+(:,)|^^^„«^■(p)e-"P^ 
d'^-p) ^ [^+{x)l^^^v^{p)e-^^'^. 

J X 

Using these one has 

r^'{x,p) = 2 (e-^^^y I?(x,y)|^^^o7V(p) , 
r''{x,p) = -2 fe-^Py D{x,y)\y^^,j%^ip} . 

J p 



(29) 
(30) 

(31) 
(32) 



is finite in the perturbative vacuum of the fermions. 

Following the existing practice in classical simulations, 
we will use a temporal discretisation step that is negligi- 
ble to the spatial lattice spacing, i.e the cut-off is three- 
dimensional, and the three-dimensional momentum inte- 
grals are implicitly regularised. We give an explicit form 
of E in the spatial Fourier space: 



E(t, k) = ~Ag' 



m? — p{p — k) 



- 1 



sin LUfft coswr -t . 

y k—p 



(37) 

We define Wp- = a/ w? + p"^ . The wave function renormal- 
isation we either get by taking the second fc-derivative at 
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scale of renormalisation, which is A; = in our calcula- 
tion, or one calculates it from the real-time behaviour 
using the formula 



6Z = Nf / dt— S(i,/c) 



k=0 



(38) 



This equation is the real-time variant of SZ = 
Nfd'^'S{ko, k)/{dko)'^ at zero momentum and makes sure 
that the coefficient of vanishes in Eq. ((36)) . 

One can perform the time integral in Eq. (|38[) under the 
assumption that oscillations of the indeterminate integral 
at large times are incoherent and they are averaged away 
when the fc-integral is carried out. One finally arrives at 



6Z=- 



(39) 



the divergence is logarithmic, as expected. An analogous 
calculation delivers the scalar mass counterterm 



dt T,{t,k) 



fc=0 



P' 



±1' (40) 



which is quadratically divergent. 

To renormalise the coupling we need to go beyond the 
linear approximation in Eq. (|35p . We renormalise the 
effective potential at zero momentum. We analyse the 
non-linear response to a static field and compensate the 
force on this static field by SV. This way we do more 
than substracting divergences. We actually alter the fi- 
nite part of the theory so that the scalar potential is 
exactly as it was before coupling to fermions. This com- 
plete renormalisation will ensure the correctness of any 
comparison with the purely bosonic classical field theory. 

A static scalar field with a Yukawa coupling is similar 
to a Dirac mass. The current J is then constant in space 
and time, but it depends on the mass M = m — g^R. 
The counterterm 5V'{^ii) based on the vacuum one-loop 
diagrams reads 



m - g^R 



3 ^/ (m - g^n) 



(41) 



Expanding this integral to linear order in $ gives the 
same counterterm as we have already found in Eq. (|40p . 
To third order in $ we find in the chiral limit for the 
coupling renormalisation 5X = 12SZ as it has been also 
derived in ^t\ . 

Of course, the integral in Eq. (|4ip would be very time 
consuming to calculate in each space-time point when 
solving Eq. (p4|) . Therefore we approximate 6V' with 
a fifteenth order polynomial fitted in the range ag^ G 
[—2.5, 2.5]. The relative precision of the fit is between 1 
and 10%, (the greatest when « 0). The fit interval is 
exceeded only by extreme excitations on coarse lattices, 
and one can extend it with little effort. 

In the rest of the paper we do not write out the R index 
for the renormalised background, and all parameters are 



understood as renormalised. For simplicity, we also hide 
the counterterms in the equations we discuss, but we keep 
them in our numerics, of course. 

Contrary to Ref. [ST] , in this approach we solve equa- 
tions with divergences, which cancel in the end result. 
This makes the final removal of the cut-off impossible, 
and such a calculation is usually error-prone close to the 
continuum limit. But in this case we solve a lattice field 
theory classically and it makes no sense to even approach 
the continuum limit. This renormalisation makes sure 
that the fermionic vacuum does not alter the bosonic 
vacuum, but the classical divergences from the closed 
bosonic loops are as dangerous as before. 

By construction, fermions have now no impact on a 
static bosonic field, but there is a damping rate for dy- 
namical fields, which is given in the real scalar case by 



dtsin{kot)^{t,k) 



(42) 



For a scalar with mass fi this evaluates for homogeneous 
mode to 



-/ifi,k) 



k=0 



2fi 



d{p/2-ujp) = 



167r/i2 



(^2_4™2)3/2 

(43) 



if /i > 2m. The damping rate is directly observable from 
the numerics. Since it is proportional to Nf, it facilitates 
the measurement of the number of doublers in a lattice 
implementation. 



III. STOCHASTIC APPROACH 



In Eq. ([35)) of the previous section a separate field has 
to be evolved for each mode p and spinor index. On a 
three-dimensional lattice with N'^ sites this means a cou- 
pled set of 4 • 4 • complex ordinary differential equa- 
tions. This relatively high price might explain the fact 
that in nearly ten years time since the equations have 
been published no calculation has been carried out be- 
yond 1+1 dimensions. 

An elegant way of performing integrals with high di- 
mensionality is to employ Monte-Carlo techniques. Im- 
portance sampling is a prominent example in statisti- 
cal field theory, though its formulation for fermionic 
fields is troublesome because of the Grassmann nature of 
these degrees of freedom. Nevertheless, there have been 
promising news to the apparently impossible simulations 
at finite chemical potential [6lj or in real time [g^I . 

In fact, the situation in our semiclassical nonequilib- 
rium setting is much simpler than in Euclidean field the- 
ory simulations. We know everything about the initial 
fermion ensemble and we will set up evolution equations 
for the members of this ensemble. At any later time an 
averaging over these members will tell the propagator 
D{x,y). 

Notice that we could formulate Eqs. ([6]l7]) as well as 
the back reaction (|12j|14p in Eq. (fTT)) without any refer- 
ence to the spectral function, which is complementary to 
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the symmetrised propagator D{x, y). We will replace the 
commutator of anticommuting operators by the product 
of plain complex numbers in D. To accept this simplifica- 
tion we have to show that the two-point function defined 
in terms of this simple product obeys the same equations 
of motion as D and that it also starts from the same 
initial condition. 

Let us introduce a set of classical spinor stochastic 
variables as c-number fields: i/'M(a;) and tpp{x). Only 
together can these "male" and "female" fields form a 
meaningful physical quantity, but the male and female 
roles are interchangeable: 

D{x,y) {iPM(x)i^F{y)) = {tpF{x)i)M{y)) ■ (44) 

The reason for why we need two spinor fields is that with 
a single spinor field only positive semidefinite correlators 
can be modelled, whereas in Eq. (I23p has negative 

eigenvalues. 

So that D in Eq. ([44|l obeys Eqs. ^ and ([7]) we require 
that both the male and female stochastic spinors follow 
the usual Dirac equation: 

{iYd^ -m + gRe^{x) - iglui^{x)j^)ipg{x) = 0. (45) 

The g (gender) index represents M or F. 

The currents expressed in terms of the stochastic fields 
read 

J^ix) = gTTD{x,x)^g{^+ixh'>ijMix)) , (46) 

jPS(x) = 5Tri?(x,a;)75=.g(^+(x)7VV'M(a;))(.47) 

Due to the interchangeability of and ipp the scalar 
and pseudoscalar currents are manifestly real and imag- 
inary, respectively. 

We have to make sure to satisfy Eq. (fT5|) . For this we 
define the Fourier transformed stochastic fields: 

Mp)- /e'f^->,(f), ^,(a;)= fe-^P^^'MP). 

J X J p 

(48) 

To reproduce Eq. (|15p we require 

{^M{mU<l)) = C^^f5{p- <1)\m{p) (49) 

To actually realise an initial ensemble with this correlator 
one has to solve the eigenvalue problem oiM.{p). This we 
have actually done already when we introduced the mode 
functions and denoted the eigenspinors as u^^-* w^^^ 
and t!^^-* corresponding to the eigenvalues +1, +1, —1 and 
— 1 respectively. 

We can express the stochastic spinor fields in terms of 
the eigenspinors as follows: 

4,mAp) = 71 E iUp)^'{p) ± ^s{p)v^{p)) (50) 

if and if are the primary complex random variables we 
use: 

(r(p)r'(p1^) = {2Trf5{p-q)5s^Al-2n%{p)), 
('7^(pV(p1^) - (2^)35(p-g)<5,,,,(l-2nl(p)). 



All other two-point correlators vanish. (Actually, these 
variables could be chosen real and do not necessarily have 
to be Gaussian.) Notice that nothing on the right hand 
side of Eq. ([50)1 bears a gender index, but the male and fe- 
male fields have different signs for the antiparticlc compo- 
nent. This allows for the stochastic representation of the 
hermitian matrix with negative eigenvalues in Eq. (j49p . 
With ^ and rj we actually simulate the ladder operators: 
this is possible since the ladder operators always appear 
in the expectation value of a commutator. 

The eigenvalues of the correlator (V'm(p)V'f(p))i which 
is a matrix in Dirac indices, actually represent the par- 
ticle number: they take the vaule ^ — n^_^\'p) for the 

fermions, and n^f\p) — ^ for the antifermions. By proper 
initialisation, one can start from a polarised fermion gas, 
or, one can set a constant non- vanishing baryon density, 
as we anticipated. In a completely symmetric setting we 
can read out the particle number by taking the deter- 
minant of the correlation matrix (in momentum space), 

which shall be (n{p) — i)*. 

At this point we return to the question of numerical 
feasibility. The expectation value in Eq. (j44p turns into 
an average over E pairs of spinor fields in practice, where 
E is finite number. The statistical error in Eq. (j44p prop- 
agates through Eqs. (|12m4p into the scalar equation. The 
statistical noise in the back reaction may induce artifical 
production of scalar fluctuations. Thus, checking for the 
_E-dependence of the final result is an essential part of 
using this scheme. If the required number of spinor pairs 
(E) turns out to be higher than the number of lattice 
sites, the standard deterministic mode function expan- 
sion is the cheaper and more precise option. This is typi- 
cally the case in 1-1-1 dimensions. Increasing the number 
of dimensions, however, one can in most cases keep E 
around the linear lattice size or less, and the stochastic 
method can be by several orders of magnitude more ef- 
ficent than the deterministic algorithm, both in memory 
need and in time. 

For future reference we give the actual form of the 
spinor equations as well as their initialisation in Ap- 
pendix 1^ Since the system we analyise is implicitly un- 
derstood to be discretised on a lattice, some comments 
on lattice doublers are due in Appendix IB] 



IV. EFFECTIVE SCALAR DYNAMICS 

In this section we present the numerical analysis of a 
real scalar field coupled to fermions as introduced above. 
For the sake of simplicity of the implementation we re- 
strict our numerics to 2-f 1 dimensions. 

We perform the renormalisation of the effective poten- 
tial as already anticipated, but no wave function renor- 
malisation is necessary. In Fig. [1] we give SV'{^) by 
evaluating the two dimensinoal variant of Eq. (I4ip on 
a large lattice for various fermion masses. In the plot we 
used a for the lattice spacing. Notice that in the massive 
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FIG. 1: Renormalisation of the effective potential. For the 
thick hne we used Wilson fermions in two spatial dimensions 
with r = 1. The other lines have been calculated in the pre- 
sense of doublers. The breaking of chiral symmetry manifests 
in the asymmetry of SV'{^) around zero. 



case with broken chiral symmetry we loose the 
symmetry. For this reason we use massless fermions and 
compensate for the doublers as detailed in Appendix IB] 

In the following we discuss a few test cases to explore 
the capabilities of this semiclassical approximation. To 
better see the effects of the fermions we always run the 
purely classical simulation with the same initial condi- 
tion (same random seed) in parallel. For reproducibility 
we give the parameters in the figure captions: the linear 
lattice size (N), the Yukawa coupling (g), the fcrmion's 
inital temperature (T/), the scalar mass and coupling 
(A), and the number of spinor fields (E) in the ensem- 
ble. In these experiments we used two-component chiral 
fermions. These parameters and the data on the plots 
are given in lattice units (with a = 1). 

In our first exercise we plot the damping of the scalar 
homogeneous mode in Fig. [2l The exponential with 
expected rate (7 = in 2-1-1 dimensions with 2- 

component spinors) nicely forms an envelope of the calcu- 
lated evolution. It was important to use a large volume, 
otherwise the damping stopped at about N/2 time and 
recurrences occur. In fact, one assumes infinite volume 
in the derivation of the decay rate. 

Let us now consider an example where the fermions 
start from a finite temperature state and transfer energy 
to the bosonic vacuum. We set up an experiment with 
a small noise in the bosonic sector, fi^ — 0.25, A = 6, 
g — 0.25 and Tf = 1. To our surprise, there was no 
boson production at all, but the small initial scalar noise 
was transformed into fermions with a rate comparable to 
7. It seems that in the semiclassical approximation the 
production of quantum fluctuation is a one-way channel 
of interaction. 

To see how the energy is transferred to fermions regard- 
less to our thermodynamical preconceptions we present 



FIG. 2; The homogeneous mode of the scalar field is exponen- 
tially damped at the expected rate. (Parameters: A'^ = 1024, 
g = 0.5, Tf = 0, p2 = 0.25, \ = and E = 20) 

the results of our third experiment. The scalar field 
is now started from a non-thermally excited state with 
an isotropic particle distribution peaked around the mo- 
menta |fco| = 0.5 with n(|A;o|) — 10. The initial energy 
density was w 1. What we see in Fig. [3] is a counter- 
intuitive anti-thermalisation, where all energy that can 
be possibly transformed to quantum fluctuations is taken 
away from the background. This also happens in the 
purely bosonic Hartree approximation, but here in the 
fermionic case the particle number is capped at a value 
close to 1/2 due to Pauli blocking. The modes above 
|fc| > 1 are quickly excited (at the order of damping 
time). The low momentum modes are filled up on a 
much slower scale. At this point we remark that a 
two-dimensional classical scalar theory comes into non- 
thermal (quasi) fixed points for a wide range of initial 
conditions. For a similar classical system we found that 
the evolution to equilibrium can be extremely slow, gov- 
erned by a power law [6^ . In this example, too, the scalar 
spectrum evolves into an approximate power law with an 
exponent of « —1.8(2). The effects of fermions manifests 
merely as an overall coefficient in the spectrum. 

We finally show an example where the classical ap- 
proximation is expected to work well. We start the clas- 
sical system from the center of a double-well potential. 
There is a rapid particle production fueled by the spin- 
odal instability. The resulting scalar spectrum is not 
far from a Bose-Einstein distribution. Of course, this 
closeness to quantum equilibrium is temporary: the slow 
classical thermalisation drives the system towards classi- 
cal equipartition. Coupling this scalar field to fermions 
switches on a dissipation, as one can see in the plotted 
energy density in Fig. |4l 

The scalar spectra in Fig. [4] are close to a quantum 
equilibrium distribution with some chemical potential. 
The UV end of the spectrum is distorted by lattice arte- 
facts, but otherwise the linear fit of \og{l /n{p) + 1) is 
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FIG. 3; Anti-thermalisation. The bosonic excitations are 
transformed to fermions until fermion-production is cut by 
Pauli blocking and a close-to-infinity temperature sets in. 
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adequate (see dotted line in Fig.d]). The closeness to the 
quantum equilibrium distribution is maintained through- 
out the time we followed the dynamics. The energy drain 
of the fermion field does not bring the scalars out of this 
equilibrium but imposes a steady cooling. 

These numerical experiments lead us to the negative 
conclusion that the energy transfer between the classical 
and quantum degrees of freedom is unidirectional. As it 
is also known, the produced quantum particles do not 
scatter on each other, and the energy transfer between 
modes through the inhomogeneous background is ineffi- 
cient. Can this semiclassical approximation be used then 
at all? 

We do think that in some circumstances this low-cost 
solution to add fermions to a classical field simulation 
is adequate. The rate of fermion production is correctly 
given account for, albeit these particles will not ther- 
malise. The very mechanism of particle production and 
the simultaneous loss of energy in the bosonic sector is 
well described as long as fermions are not created in such 
an abundance that their non-thermal distribution could 
have impact on the back-reaction. In fact, bulk observ- 
ables, such as the scalar effective potential, prethermalise 
[63 |. i.e their value before thermalisation can be used as 
an estimate to what one would find after equilibration. 
Even if the fermion distribution is non-physical, the evo- 
lution of the bosonic background can be well approxi- 
mated. If, however, the back-scattering of the produced 
fermions to bosons becomes relevant, this semiclassical 
approximation will no longer be applicable. One can ac- 
tually monitor the fermion particle numbers to check for 
relevance of (the absence of) back-scattering. 
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FIG. 4: Scalar field with a spinodal instability. Top: en- 
ergy density of the scalar field with or without coupling to 
fermions. For comparison, three different Bottom: the par- 
ticle spectra at t = 150. At and around the time shown, 
the scalar spectrum is close the Bose-Einstein distribution, 
especially when coupled to fermions. The dotted line in the 
inset plot is the thermal fit (/3 = 4.9). (Parameters: N = 64, 
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V. FERMIONIC DECAY OF OSCILLONS 

We can consider the semiclassical approximation safe 
if the resulting fermion energy density is small. However, 
one of the justifications for the classical approximation is 
the high bosonic occupancy, which will inevitably gener- 
ate an energy transfer into the fermionic fields. In such 
cases the approximation will break down within a short 
time, which is likely to be the damping time we discussed 
in the previous section. 

In the this section we turn to applications where clas- 
sicality has an other justification. If the particle content 
is very low and but there is dilute network of classical 
structures, such as topologial defects, their evolution can 
be well described by the non-linear wave equations. The 
decay of these structures into particles is mapped to the 
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production of classical waves ( "ripples" ) by the classical 
equations. As this mostly happens in the ultraviolet, the 
classical approach is not justified for describing particle 
production, in contrast to its usefulness in the case of the 
macroscopic networks, like cosmic strings. 

We addressed this deficiency of the classical approx- 
imation in Ref. [65|, where we introduced a stochastic 
approach to the bosonic mean-field approximation, sim- 
ilar to the method presented in this paper. We found 
that mimicking the quantum distribution by an anal- 
ogous classical noise (following the so called "just-the- 
half" prescription) introduces undesired time-dependent 
renormalisation effects to the effective potential. Instead 
we solved the inhomogeneous mean-field equations and 
found that on the macroscopic level, the decay channel 
into quantum particles plays negligible role, whereas it 
on microscopic scale we found deviations. Our numer- 
ical analysis suggested that oscillons, which are one of 
the classical decay products of topological defects [6^ . 
are the primary sources of quantum particles, while the 
direct radiative decay of a defect network is suppressed 
as predicted in Ref. [63]. 

To better understand how oscillons decay quantum 
mechanically we reproduce one of the experiments in 
Ref. [g^ , but we also add fermions. In two dimensions os- 
cillons are particularly stable (68l.[69j localised structures, 
when several oscillons are created in volume, they behave 
as molecules in a gas. When oscillons collide, the coalesce 
with some probability. Being this mechanism their only 
decay process (in 2-fl dimensions), the number density 
of oscillons obeys the equation fi{t) ^ n^{t). Thus, the 
classical solution is n{t) ^ 1/t, which is approximately 
manifest in classical simulations ^66j . 

We put 16 incoherent oscillons with small random ve- 
locities in a box with N = 128. We estimate the number 
of oscillons by counting the sites with an energy density 
beyond a thresold (e > 0). This number we normalise to 
its initial value and plot in Fig. [51 It takes long before the 
expected power-law solution sets in (and even then finite 
volume effects can distort it). But a small coupling to the 
fermionic fields introduces a new time scale, and the slow 
classical behaviour is replaced by a close-to-exponential 
decay. (The oscillon damping rate is about four times 
stronger than for the homogeneous mode.) This process 
reduces the amplitude of most oscillons below the thresh- 
old. After t > 100, however, the plotted estimate can be 
best fit by a power law with an exponent of —2. For the 
semiclassical evolution of these localised objects a sur- 
prisingly small spinor ensemble already provides results 
that are insensitive to an increase in E. 

In the inset plot we estimated the oscillon decay rate by 
the inverse time necessary to radiate away 100exp(— 1) 
percent of the oscillons. The rate cuts off at about ^cut ~ 
0.45. One can explain this by simple kinematics. The 
effective mass of the produced fermions nif ~ gv, where 
V = 6/i2/A is the vev of the background. The scalar 
mass in the broken phase is rrif, = —2/1^. It is not 
this bosonic mass that enters the kinematical relation 
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FIG. 5: Estimated "molecule" number in a gas of oscillons. 
When the scalar background is coupled to fermions, the slow 
classical evolution is replaced by an approximately exponen- 
tial decay. An estimate of the rate as a function of the Yukawa 
coupling is shown in the inset plot. The obtain the same 
curves using the mode function expansion would have required 
three orders of magnitude more computational resources. For 
g = 1/4 we explicitly check for the insensitivity to doubling 
the ensemble. (Parameters: iV = 128, T/ = 0, = -0.25, 
A = 3 and i? = 10 or 20. For each coupling we averaged 30 
runs. ) 



but the oscillon frequency uJosc, so the condition for the 
decay is ^Wosc > "ly. From ^cut we can tell the oscillon 
frequency: ujosc/mb = 2(?cut\/3/A « 0.9. This estimate 
is in harmony with direct measurements [66j . 

In this paper our aim is not to explore the parameter 
space, and to analyse the mechanisms oscillon decay. In- 
stead, we put forward a low-cost technique to check exist- 
ing and future analyses of defect evolution for fermionic 
quantum corrections, complementing work already done 
for bosonic ones [11]. We plan to investigate the evolution 
of cosmic strings for such contributions from quantum 
degrees of freedom in a future publication. 



VI. CONCLUSIONS 

In this paper we propose a low-cost integration scheme 
for the fermionic path integral, which leads to equations 
that are equivalent to the mean-field approximation stud- 
ied earlier by Aarts and Smit. These equations also follow 
from the large- A^/ expansion of the 2PI effecitve action. 
The computational efficiency of this scheme allowed us to 
do simulations beyond l-fl dimensions. This stochastic 
method is a generalisation of our earlier technique devel- 
oped for scalars in Ref. [6^ . 

We calculate several test cases on a scalar example and 
study what the irreversible phenomena can be captured 
by this simple method. We confirm that the fermions, 
once created, can no longer scatter on each other, but 
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this is not the only obstacle that hinders thermalisation. 
The damping of the classical oscillations are correctly 
given account for, but the back-scattering of the fermions 
into bosons is absent. In the language of the mode func- 
tion expansion, fermions (and also other quantum fluctu- 
ations on the Hartree level) are represented by far more 
dynamical variables that the background. If these vari- 
ables strive for classical equipartition (as usual in a cou- 
pled set of non-linear differential equations), the energy 
left in the background is negligible. This suppression of 
the background becomes stronger with higher dimension- 
ality, and was less relevant in former 1+1 dimensional 
calculations. 

As it was remarked in Ref. [Slj the fermion spectrum 
can become close to thermal, and this raised hope that 
the inhomogeneous Hartree approximation is still capa- 
ble to account for an approximate thermalisation. It is, 
however, more likely, that it is the praticle production 
mechanism that brings the fermions close to equilibrium, 
rather than scattering. 

Even though scattering cannot drive the fermions to- 
wards equilibrium, we expect that the back-reaction of 
the often not-far-from-thermal fermion field has an ap- 
proximately thermal back-reaction due to prethermalisa- 
tion of the fermionic current [64] , and the lack of thermal- 
isation has little impact on the background field. This 
assumption becomes even more plausible if we assume 
that the fermions leave the scene, once created. 

In some physical situations it is difficult for a fermion 
to leave. If we apply the presented scheme to the Yang- 
Mills equations, and solve the semiclassical chromody- 
namics, it will be difficult for fermions to be reabsorbed 
by the plasma. This puts jet-quenching outside of the 
range of validity. But a semiclassical simulation of the 
freeze-out of the plasma is not ruled out by the aformen- 
tioned deficiencies. 

The numerical calculation of the fermion spectrum in 
baryogenesis scenarios is a more viable application. If 
baryogenesis is driven by a first order phase transition, 
the presented equations can give account for CP violation 
as well as the departure from equilibrium without relying 
on gradient expansion, and thus, allowing for thin walls. 
For the subsequent thermalisation, however, one has to 
make further assumptions. 

The scheme is best applicable for systems with low par- 
ticle numbers and genuine inhomogeneities, like a dilute 
network of topological defects, such as cosmic strings. In 
this context fermion production is local, and the pro- 
duced particles spred in space. This results in small 
praticle numbers, and we expect that the lack of ther- 
malisation will introduce very little distortion into the 
back-reaction. 

In conclusion, for cases where the inhomogeneities in 
the background are ment to be "particles" , a big vol- 
ume is less relevant, and other techniques, such as the 
2PI effective action on a homogeneous ensemble may be 
more favourable. For the large-scale classical simulations 
with inhomogeneous classical structures, however, the in- 



homogeneous 2PI approach would be beyond feasibility. 
In such situations, the Hartree approximation already 
includes the leading quantum corrections, as well as en- 
dorses fermions. The technique in this paper has made 
these type of calculations affordable. 
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APPENDIX A: REPRESENTATION OF THE 
SPINOR FIELDS 

The initial conditions for the spinor fields is given 
in terms of the u''{p) and f^(p) eigenspinors. Here we 
present the actual form of these eigenvectors. We use the 
normalisation factors for a theory discretised in a volume 
V. 

The formulas are based on a naive fermion action. 
They can, however, be easily rewritten for the Wilson 
fermions by replacing to to m -f ip^ when discussing the 



p mode, with pj = j^smpja/2. 



Chiral basis in 3+1 dimensions 



In the chiral base we define the gamma matrices as 



1 

1 



Consequently, 



r = 



a' 
-cr* 



a' 

-a' 



7 = 



1 

-1 



-1 
1 



The hermitian 4-by-4 matrix in Eq. 
basis reads 



V 



pa TO 



(Al) 
in the chiral 

(A2) 



Here pij stands for the combination of the Pauli ma- 
trices: 



pa 



P3 Pi - tP2 
Pi + ip2 ~P3 



(A3) 



We used the standard notation pj — asinpja. 
In the chiral base the eigenvectors are given as 



1^ +^p\^ 




Xp u 



(2) 



\P\ 



P \ —TO 



8 Xp 

^Xp 
(A4) 
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with = 2^^(J"-+|p1) ■ denote the eigenvectors of 

Picr' for the eigenvalues |p| and — respectively. 
Let us now diagonalise pjcr^: 



One finds that the eigenvectors are 
Q* 



P3 + \P\ 

Pi + ip2 
Pi - iP2 

\P\ -P3 



Xp = P+ 
Xp = P- 



-Pi + ip2 

P3 + \i>\ 

P3 ~ \P\ 

Pi + ip2 



if P3 > 
if P3 < 



u{p) = (3 



for p2 > 0, and 



uip) = P 



v{p) = /3 



v{p) = (} 



(A9) 



(AlO) 



with = l/2|p|(|p| ips). The two cases we handle 
separately for numerical stability (e.g. to avoid divisions 
by zero). 

We actually solve the Dirac equation ([IS]) for the '4}g 
field instances. This reads in chiral base 



/ ds di-id2 -iM* 
di + id2 -83 -iM* 

-iM^ -^3 -81+182 

\ -^Af, -di-id2 83 

(A5) 

with Mx — m — g^{x). We can use this equation for both 
ipM and iJ^F- For Wilson fermions ~ ^A. 



2. Majorana basis in 2+1 dimensions 

Here we work out the implementation details for a 2+1 
dimensional setting with a real scalar background. In 
2+1 dimensions we have 2-by-2 gamma matrices. In Ma- 
jorana basis all these are imaginary: 



7 



-i 

1 



7 



i 
-i 



7 



i 

1 



In 2+1 dimensions one of the gamma matries can be 
easily expressed by others: 



7*^7"'" = 47^, 7^7^ = 



(A6) 



The Dirac equation in this basis (with = m — g^) 
reads: 



-82 81 - 
81 + Af, 82 



(A7) 



The advantage of the Majorana basis becomes apparent 
with the form of this equation: the spinor field equa- 
tion is real. Although the spinor fields themselves are 
complex, their real and imaginary part follow a separate 
equation of motion. This facilitates numerical optimisa- 
tions, such as vectorised arithmetics, and it requires a 
smaller memory-to-cache bandwidth. 

We have to initialise the spinors in terms of eigen- 
spinors. For this, we diagonalise the vacuum correlation 
matrix 



Mip) 



1 



-P2 Pi - im 
pi + im p2 



(A8) 



with Q ^ pi — im, s = p2 + LUp and (3^^ — \Q\'^ + 



APPENDIX B: FERMION DOUBLING 
PROBLEM IN THE SEMICLASSICAL THEORY 

The problem of fermion doubling inevitably arises in 
any lattice implementation. Since almost all numerical 
^^jialyses of classical field theories use lattice discretisa- 
tion, an extension that incorporate fermions will also 
share this heritage. Time discretisation, however, is not 
an intrinsic parameter of the classical theory. Whereas 
the space- like continuum limit simply does not exist, we 
can always assume that our equations are in the time- like 
continuum limit. Indeed, the time-step (at) in our numer- 
ics was much smaller than the lattice spacing a = 20at. 

There are several remedies in the literature for the 
problem of doublers. We made a version of our numer- 
ics using Wilson fermions, but the explicit breaking of 
chiral symmetry introduced a linear term in the poten- 
tial. Although this can be renormalised away, not only 
the vacuum, but also the physical excitations will also 
contribute and introduce artefacts in the scalar effective 
potential. This effect will vanish in the continuum limit, 
but in a semiclassical theory, we cannot go close to the 
continuum limit, by construction. 

The other low-cost solution could be the use of stag- 
gered fermions. These are, however, special to two or 
four dimensions, and some of the doublers will be kept. 
To avoid complications on the level of the equation of 
motion we dropped this idea too. 

In the presented numerics we simply used the naive 
fermion discretisation and introduce an effective fiavour 
number, in which we compensate for a pair of doublers 
in each spatial direction. We could do this since in our 
simple model there are no anomalous diagrams where 
doubling fermions could cancel. 

There is, however, a time-like discretisation, too, which 
can be a source of time-like doublers. To eliminate them, 
Aarts and Smit used a linear combination of two different 
fiavours, and the two degrees of freedom have both been 
made physical. 

In the following we analyse the real-time Dirac- 
equation to understand how such doublers affect our nu- 
merics. 

The free Dirac propagator on spatial lattice in momen- 
tum space reads 



Dit,p) 



m +Pjj' 



7 

cos{ujt) — i-^ siniiut) 



(Bl) 
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with pj = sm{apj) and w?. = + J^jP]- 
time-like continuum limit lo = up must be satisfied so 
that Eq. (|Bip solves the Dirac equation. If time is discre- 
tised as the average of the forward and backward deriva- 
tive, then Dirac equation takes the following form: 



m)D{t,p) = 0. 



(B2) 



If we use the Aat in the equations, the Xe() func- 
tion will always give one in the source J, since there 
we close the Fermion loop by evaluating the propagator 
equal space and time. At that point we need to compen- 
sate for the extra factor two, compared to the continuum 
limit. We achieve this by removing a factor of two in 
Eq. (|B5P from the initial value of D. 



Inserting Eq. (|B1|1 into Eq. (|B2p we get following con- 

sin(a;at). For an extremely 



straint: u; 



top with CO ■■ 



anisotropic lattice {at «C a) either w « w or cD « Tr/a-t — uj, 
since up is limited by the spatial cut-off. This means, 
that there are two solutions (the doublers) which can be 
worked out explicitly as 



Di{t,p) 

D2{t,P) 



2ujp 



7 

{ujpt) — i— sin{ujj;t) , (B3) 



-PjT 



2ujp 



COii{uJpt){-lY 



(B4) 



where s is the index of the time-slice i, i.e. t — ats. The 
sum of these solutions is the standard lattice propagator: 



Aat(i,p) = 2 



- PjT 



cos(wjji)xe(s) - i— s\n{ujpt)xo(s) 



(B5) 

where we introduced the Xe() and Xo() functions, which 
is one if their integer argument is even or odd, respec- 
tively, and zero otherwise. Indeed, Fourier transforming 
Eq. ((B5)) yields (in the ot/a < 1 limit) 

Aat(p) = 7r5(po - LUp) [(to +pj7-') + 7°tjp-sgn(po)] • 

(B6) 

We get the continuum propagator from Eq. (jB6[) by re- 
moving the bars. The staggered nature of the lattice 
propagator is also manifest in spatial coordinates: e.g. 
TrZ?(i, af)7^ is only then non- vanishing if xi/a is odd. 



At zero time we start our system with excitations de- 
scribed by the Di propagator. The space and time- 
dependence of the background will result an inhomo- 
geneous propagator Di{x,y). Had we started from an 
initial condition corresponding to the D2 propagator, 
the evolution would have lead to D2{x,y). Inserting 
Aat = Di{x,y) +D2{x,y) or Aat = Di(x,y) - D2{x,y) 
into the inhomogeneous Dirac equation one discovers that 
these linear combinations decouple: they communicate 
only through the back-reaction to the scalars. The vari- 
ous Lorentz-components of D\^i{x, y) and l)iat(a;, y) cou- 
ple to the background at different time slices, depending 
on the parity of x'^ — y'^ . If the background is a smooth 
function of time, Diat and Z)iat will evolve on the same 
background, up to an error at. Thus, their differ- 
ence, D2{x,y) is suppressed by the time- like spacing, i.e. 
if there is no D2 component in our initial condition, the 
production of doubler particles will be small compared to 
t^e standard particles. In the back-reaction and the mea- 
sured spectra both types of excitations contribute indis- 
tinguishably. (Notice that Di{t,p) and D2{t,p) are iden- 
tical at equal time, where these observables are taken.) 
Similar ideas have been implemented to tackle the species 
doubling problem in the context of the hard thermal loop 
effective action of the electroweak theory in Ref. [t^ . 

To check these ideas we plotted the damping of the 
scalar field in Fig. [51 For this calculation in 2-1-1 dimen- 
sions we used the effcctvc flavours number 1/4. An er- 
ronous estimation of the number of flavours should have 
generated an unexpected factor 2 in the rate. 
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